2.7 ms
83.9 ms
565 s

Statistical inference with Julia

Group Retreat Heidelberg, Matthias and Maurice, 15. September 2020


395 ms

Poisson parameter:

9.8 ms
254 ms

(Package imports above title)

24.3 μs

Before we start

Tutorial requirements (the complete process should take <30 min.)

  1. Julia installation

  2. Optional: Atom and Juno (the Julia editor)

  3. Julia packages listed below:

1.6 ms
Status `/private/var/folders/29/l75f_tqn5fzb5qylqt0w2xbr00012v/T/jl_q64l4i/Project.toml`
  [6e4b80f9] BenchmarkTools v0.5.0
  [0c46a032] DifferentialEquations v6.15.0
  [31c24e10] Distributions v0.23.12
  [4138dd39] JLD v0.10.0
  [91a5bcdd] Plots v1.6.6
  [c3e4b0f8] Pluto v0.11.14
  [7f904dfe] PlutoUI v0.6.4
  [f3b207a7] StatsPlots v0.14.13
  [fce5fe82] Turing v0.14.4
367 ms

How to install and precompile Julia packages

  1. Start Julia in the terminal (Applications -> Julia-1.5)

  2. Enter “]” to get into Julia package mode (pkg>)

  3. Install a package X by writing “add X” (e.g., "add Plots", repeat for each package)

  4. View all installed packages with “status”

  5. Precompile all packages with “precompile”

  6. Leave the package mode by pressing backspace (julia>)

  7. (General info: Import a package by “using X” for a package X)

How to start this Pluto notebook

  1. Start Julia in the terminal (Applications -> Julia-1.5)

  2. Import Pluto by typing "using Pluto"

  3. Open Pluto in a browser by typing "Pluto.run(1234)"

  4. Open this notebook file via the Pluto user interface

(Currently, Pluto is supposed to work best with Firefox or Chrome)

10.9 μs

Content of this Tutorial

  1. Introduction to Julia

  2. Performance tips and a Gillespie algorithm

  3. Differential Equations in Julia

  4. Bayesian Inference and Probabilistic Programming

  5. Everything from above, just combined :)

  6. Further Information/References

5.3 μs

1. Introduction: What is Julia?

Goal: Write code as in Matlab/Python/R which is as fast as C/C++/Fortran

Basic language principles

  • Multi-platform via LLVM

  • Multi-paradigm (functional, object-oriented, ...)

  • Multiple dispatch

  • Dynamic with type inference

  • Just-in-time (JIT) compilation

  • Broad metaprogramming features

  • Open source

10 μs

The Two Language Problem

Effort and language barriers

  • prototype code with a language that is fast to code with (Python, R, Matlab)

  • tranfer the prototype code to a language that runs fast (C, C++, Fortran)


The numpy package in Python (basic array class for fast computing)



A typical Julia package (Turing.jl)


Julia solves the two language problem

  • you can prototype code in Julia

  • you can improve code performance in the same language

  • no language barrier

25.7 ms

Basic syntax

Getting help

Precede unknown object/expression with a question mark (?) to see the docs (Pluto has Live docs for this)

5.9 μs
Some basic data structures
  1. Tuples (ordered, immutable)

  2. Named Tuples (ordered, immutable)

  3. Arrays (ordered, mutable)

  4. Dictonaries (unordered, mutable)

Tuples and Named Tuples are constructed using paratheses

9.5 μs
3.1 μs
4 ms

Elements are accessed by position (In Julia indexing starts at 1!)

6.6 μs
3.9 μs
2.8 ms

Named Tuples can additionally be accessed by name

7.3 μs
3.1 ms

To form collections of related values, use Dictionaries or Arrays

6 μs
4.5 μs
40.9 ms

Arrays are accessed by indexing

6 μs
2.4 μs

Dictionaries are (unordered) collections of (key, value)-pairs, values are retrieved via corresponding key

5.7 μs
12.7 ms
Loops and control-flow
2.9 μs
2.5 μs
7.5 ms
2.9 ms
595 ns
Functions
2.9 μs
27.8 μs
11.2 ms

Type inference

2.8 μs

Julia is a dynamic language, you don't have to specify types, but you can. Either way, Julia will try to infer types that allow fast subsequent computing.

2.6 μs
some_array
25.1 ms
Array{Any,1}
5.5 μs
744 ns
11.4 ms
2.1 μs

Abstract and concrete types

4.3 μs

Relationships between types

10.6 μs
10.1 μs
6 μs
5.4 μs
16.3 μs

Parents and children of types

7.1 μs
4.5 μs
117 ms
52.1 ms
46 ms

Build your own types

9.4 μs
5.1 μs
3.1 ms
3.7 μs
2.4 ms
1.3 μs
591 ns

Multiple Dispatch

6.1 μs

Multiple dispatch is the primary paradigm in Julia; it allows to define the same function for multiple type combinations for the input variables. The most appropriate (i.e. fastest) version is then dispatched to a certain input.

8 μs

+ is a function in Julia with many different versions (called methods)

4.9 μs
4.2 μs
4.1 μs
65.4 μs
9.6 ms
10.7 ms

Create the function mysum with two different methods

7.2 μs
37.2 μs
506 ns
7.9 ms
6.2 ms

2. Performance tips

3.5 μs
For-loops vs. loop-fusion/vectorisation/broadcasting
  1. Exercise: Which version is the fastest?

3.2 μs
2.9 ms
47.8 ms

Check box to show solution

70.7 μs

You will do it!

5.2 ms
In-place functions
  1. Exercise: What could 'in-place' mean? Which version is in-place and what might be faster?

3 μs
28.1 μs

Check box to show solution

60.2 μs

You will do it!

8.9 μs
Type-stable code
  1. Exercise: What could 'typle-stable' code mean? Which of the following versions has a type-instability? Which is faster?

5.3 μs
37.9 μs

Check box to show solution

59.1 μs

You will do it!

8.5 μs

More performance tips:

  • Use macros and packages to help you find bottlenecks (not guess about runtime, just test it!); such as @btime, @allocated and @profview

  • Use the @code_warntype macro to find type-instabilities

  • Copying arrays can be expensive, often Views of arrays are enough and much faster!

  • Access multi-dimensional arrays along columns since Julia stores in column-major order

More resources:

4.5 μs

2. Gillespie algorithm

Simple division process, Julia says 'Challenge accepted!'

In the following we have a stochastic simulation algorithm for a simple division process of a single cell with exponential division times (Gillespie algorithm).

Let's check out how Julia code looks like and how fast it is!

751 μs
35.2 μs

Let's check if we get expected results! We have a linear Gillespie algorithm, thus we know the analytical result for the mean cell numbers as an ordinary differential equation: dxdt=λxx(t)=x0exp(λt)

2.6 μs
80.4 ms
4.2 ms
39.1 ms
10.6 s
Exercise (optional): Transfer the code to an editor of your choice (Python, R, Matlab) and compare the runtimes with the Julia version!
3.1 μs

Check box to show solution

70 μs

You will do it!

709 μs

3. DifferentialEquations.jl

5.1 μs

An ordinary differential equation (ODE) model

Below you find an ODE model for the famous Lotka-Volterra system; an example of pre-ancient theoretical biology (starting around 1910-ish). It describes predator-prey relationships.

1.4 ms
29.9 μs

α =

β =

γ =

δ =

125 μs
12.1 ms
1.4 s
789 ms
Exercise:
  • Which variable is prey? Which is predator?

  • What do the parameters α, β, γ and δ mean biologically?

4 μs

Check box to show solution

70.6 μs

You will do it!

9.3 μs

What else does DifferentialEquations.jl offer? Basically everything...

  • Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)

  • Ordinary differential equations (ODEs)

  • Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)

  • Stochastic ordinary differential equations (SODEs or SDEs)

  • Random differential equations (RODEs or RDEs)

  • Algebraic differential equations (DAEs)

  • Delay differential equations (DDEs)

  • Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)

  • (Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)

And a whole ecosystem around it (SciML):

  • Parameter estimation

  • Sensitivity analysis

  • Bayesian inference

  • Neural networks + Differential equations

  • ...

8.5 μs

4. Turing.jl

2.8 μs
Bayesian inference with probabilistic programming (Poisson example)

Let's load some data: (two data sets available data5 and data50)

1.6 ms
2.1 s
3.5 μs
1.7 s

Define a prior:

3.2 μs
7.3 ms
637 ms

Set up the Turing model:

2.8 μs
58.9 μs

Run the inference:

3 μs
11.8 s
956 ms
1.3 s
1 s
Exercise: Change the Gamma prior to an Uniform prior! Which kind of uniform prior is a good choice?
2.6 μs

Check box to show solution

71.1 μs

You will do it!

9.4 μs
Some theoretical background for this example

Using Bayesian inference: How to learn something about a parameter θ of a model M given some data D? We first interpret θ and D as random variables. We then incorporate (subjective) prior knowledge about θ as a probability distribution p(θ) (the 'prior'). Our goal is now to calculate the 'posterior' p(θ|D) (the new information for θ given the data D). To get this we use Bayes' rule: p(θ|D)=p(D|θ)p(θ)p(D) (in words, posterior = likelihood x prior / evidence ).

What is probabilistic programming? Probabilistic programming is a programming paradigm in which probabilistic models are specified and inference for these models is performed automatically. It represents an attempt to unify probabilistic modeling and traditional general purpose programming in order to make the former easier and more widely applicable (ref: Wiki).

Constitutive mRNA expression at steady state: The stationary/steady state distribution for mRNA numbers X with linear, constant synthesis and degradation rates γ and δ can be shown to be the Poisson distribution XPoi(λ) with λ=γδ. Further, to get an analytical result for the posterior p(λ|D), we will make use of the Gamma distribution as its a conjugate prior to the Poisson likelihood; we have p(λ|D) (Gamma) ∝ p(D|λ) (Poisson) x p(λ) (Gamma).

27.3 μs

5. Turing.jl and DifferentialEquations.jl

5.5 μs
Simple division process, continued


What we have learned so far:

  • Principles of Julia, basic syntax, performance tips

  • How to Gillespie-simulate a simple cell division process

  • How to run Differential Equations in Julia

  • How to estimate model parameters using Turing

Now let's bring it all together! Goal of this last part:

  • Use our Gillespie algorithm to create some in silico data for the division process

  • Use DifferentialEquations.jl to build an ODE model for the mean cell numbers

  • Use Turing.jl on top of that to estimate our ODE model parameters, given our Gillespie data!

893 μs

1. Create some in silico data (with our 'true' parameter value for λ):

7.1 μs
351 μs
114 ms

2. Build an ODE model:

Exercise: Complete the function simple_div_ode (below) with the correct ODE model!
16 μs

Check box to show solution

75 μs

You will do it!

8.2 μs
19 μs
2.3 ms
58.8 ms

3. Setup the Turing model:

Exercise (for the experts): Complete the Turing model below, assuming a Normal/Gaussian error model for our likelihood!
9.1 μs

Check box to show solution

73.5 μs

You will do it!

9 μs
11.9 ms
1.3 s
1.1 ms
105 μs

4. Bayesian inference:

Enjoy fast Julia code for fitting an ODE model with full Bayesian inference in just under one second! :)

(Note: We use a gradient-based sampler; automatic differentiation is carried through the complete Turing model including our ODE, which makes this fast)

15.7 μs
12.2 s
3.3 ms
3.2 ms
240 ms
440 ns

6. Further Information and References

Some valuable links


Editors/Programming Environments for Julia

Editors

Notebooks


Packages and Ecosystems

Larger ecosystems

  • SciML - Open Source Scientific Machine Learning

  • The Turing Language - Bayesian inference with probabilistic programming

  • Flux - Pure-Julia approach to machine learning with Julia's native GPU and AD support

  • JuliaStats - Statistics and Machine Learning made easy in Julia

  • BioJulia - Bioinformatics and Computational Biology in Julia

  • JuliaGraphs - Graph modeling and analysis packages for Julia

  • JuMP - Modeling language for Mathematical Optimization

  • JuliaIO - Group for a unified IO infrastructure

More particular packages (sometimes part of ecosystems above)

  • DifferentialEquations - Multi-language suite for high-performance solvers of differential equations

  • DataFrames - In-memory tabular data in Julia

  • Distributions - Julia package for probability distributions and associated functions

  • Zygote - next-gen automatic differentiation (AD) system for source-to-source AD in Julia

  • KissABC - Pure julia implementation for efficient Approximate Bayesian Computation

  • NestedSamplers - Implementations of single and multi-ellipsoid nested sampling

Performance packages

Notebooks/Visualisation

  • Plots - Powerful convenience for Julia visualizations and data analysis

  • Pluto - Lightweight reactive notebooks for Julia

  • IJulia - Julia kernel for Jupyter

17.9 μs




9.5 μs

Well done, that's all!









10.3 μs

Presenter's notes

Time management

  1. Introduction to Julia 18 min.

  2. Performance tips and a Gillespie algorithm 18 min.

  3. Differential Equations in Julia 18 min.

  4. Bayesian Inference and Probabilistic Programming 18 min.

  5. Everything from above combined :) 18 min.

  6. Further Information/References 0 min.

  7. Total = 5 * 18 min. = 90 min.

Example formula dxdt=αx

Example table

MeaningVariable
Number of peoplepeople
Average number of slices each person eatsavg
Number of slices on a piece of pizzaslices

Example box

14 μs

Solution

Yes that is right, that's a lot of pizza! Excellent, you figured out we need to round up the number of pizzas!

132 ms

Example input types

a =

b =

c =

d =

e =

f =

172 μs
6.5 μs